Install and load required ackages

library(clusterProfiler) #BiocManager::install("clusterProfiler")
library(org.Mm.eg.db) #BiocManager::install("org.Mm.eg.db") ## this is for mouse! 
library(DOSE) #BiocManager::install("DOSE")
library(ggnewscale) #install.packages("ggnewscale")
library(pathview) #BiocManager::install("pathview")
#library(ggupset) #install.packages("ggupset")
library(enrichplot) #BiocManager::install("enrichplot")
library(knitr)

#remotes::install_github("YuLab-SMU/clusterProfiler") 
#library(BiocManager)
#BiocManager::valid()

Functional analysis

The goal of functional analysis is provide biological insight, so it’s necessary to analyze our results in the context of our experimental hypothesis: Depletion of TCF7 results in an enrichment in genes required for myeloid cell identity. Therefore, we may expect an enrichment of processes/pathways related to myeloid/innate immune responses.

GO Ontologies

The Gene Ontology (GO) knowledgebase is the world’s largest source of information on the function of genes. GO terms are organized into three independent controlled vocabularies (ontologies) in a species-independent manner:

  • Biological process: refers to the biological role involving the gene or gene product, and could include “transcription”, “signal transduction”, and “apoptosis”. A biological process generally involves a chemical or physical change of the starting material or input.

  • Molecular function: represents the biochemical activity of the gene product, such activities could include “ligand”, “GTPase”, and “transporter”.

  • Cellular component: refers to the location in the cell of the gene product. Cellular components could include “nucleus”, “lysosome”, and “plasma membrane”.

Each GO term has a term name (e.g. DNA repair) and a unique term accession number (GO:0005125), and a single gene product can be associated with many GO terms.

clusterProfiler

We will be using clusterProfiler to perform over-representation analysis on GO terms associated with our list of significant genes.

Over-representation (or enrichment) analysis is a statistical method that determines whether genes from pre-defined sets (ex: those beloging to a specific GO term or KEGG pathway) are present more than would be expected (over-represented) in a subset of your data. In this case, the subset is your set of under or over expressed genes.

Clusterprofiler will take as input a significant gene list and a background gene list and performs statistical enrichment analysis using hypergeometric testing. The basic arguments allow the user to select the appropriate organism and GO ontology (BP, CC, MF) to test.


Prepping the data

resdata <- read.csv("TCF1KOvsWT_all_matrix_symbol.csv", header = T)
head(resdata)
##             input_id     baseMean log2FoldChange      lfcSE        stat
## 1 ENSMUSG00000000001  3904.959106     0.07596003 0.07628721  0.99571117
## 2 ENSMUSG00000000028    66.299713     0.30166718 0.32777276  0.92035463
## 3 ENSMUSG00000000037     4.938477     0.10152331 1.20464180  0.08427676
## 4 ENSMUSG00000000056   833.776305    -0.22561867 0.11117532 -2.02939519
## 5 ENSMUSG00000000078 25622.396565     0.35991087 0.21278531  1.69142726
## 6 ENSMUSG00000000085   377.598537    -0.14503206 0.15036558 -0.96452961
##       pvalue      padj  WT_CD8_rep1  WT_CD8_rep2  WT_CD8_rep3 TCF1_KO_CD8_rep1
## 1 0.31939050 0.5934865 3.601914e+03  4013.918639  3795.816820      4075.603038
## 2 0.35738747 0.6294404 5.894950e+01    53.379113    64.922751        78.882639
## 3 0.93283639        NA 9.991441e-01     9.531984     5.410229         4.533485
## 4 0.04241805 0.1492914 9.471886e+02   850.253010   899.180096       737.144665
## 5 0.09075523 0.2604725 1.800158e+04 25351.265763 24233.498716     32210.411110
## 6 0.33478047 0.6091535 3.866688e+02   381.279377   420.915833       324.597528
##   TCF1_KO_CD8_rep2 TCF1_KO_CD8_rep3 symbol
## 1      3943.778610      3998.723034  Gnai3
## 2        56.672210        84.992064  Cdc45
## 3         5.060019         4.096003  Scml2
## 4       767.098842       801.792608   Narf
## 5     27000.260027     26937.364441   Klf6
## 6       388.609440       363.520276  Scmh1

Create Gene lists

To perform the over-representation analysis, we need a list of background genes and a list of significant genes. For our background dataset we will use all genes in our resdata results table. For the significant gene list we will use genes with p-adjusted values less than 0.05 and absolute log2FC of 1.

first lets prep background gene set

# we want the log2 fold change
original_gene_list <- resdata$log2FoldChange
head(original_gene_list)
## [1]  0.07596003  0.30166718  0.10152331 -0.22561867  0.35991087 -0.14503206
# name the vector
names(original_gene_list) <- resdata$input_id
head(original_gene_list)
## ENSMUSG00000000001 ENSMUSG00000000028 ENSMUSG00000000037 ENSMUSG00000000056 
##         0.07596003         0.30166718         0.10152331        -0.22561867 
## ENSMUSG00000000078 ENSMUSG00000000085 
##         0.35991087        -0.14503206
#original_gene_list
# omit any NA values
background_gene_list<-na.omit(original_gene_list)
# sort the list in decreasing order (required for clusterProfiler)
background_gene_list = sort(background_gene_list, decreasing = TRUE)
head(background_gene_list)
## ENSMUSG00000100426 ENSMUSG00000079481 ENSMUSG00000094840 ENSMUSG00000018925 
##          11.002034           8.518602           8.198348           7.755458 
## ENSMUSG00000091387 ENSMUSG00000023132 
##           7.649971           7.247676

now lets prep significant gene set

# Extract significant results (padj < 0.05) using subset()
sig_genes_df = subset(resdata, padj < 0.05)

# From significant results, we want to filter on log2fold change
sig_genes <- sig_genes_df$log2FoldChange

# Name the vector 
names(sig_genes) <- sig_genes_df$input_id

# omit NA values
sig_genes <- na.omit(sig_genes)

# filter on min log2fold change (log2FoldChange > 1)
sig_genes <- names(sig_genes)[abs(sig_genes) > 1]
head(sig_genes)
## [1] "ENSMUSG00000000301" "ENSMUSG00000000303" "ENSMUSG00000000489"
## [4] "ENSMUSG00000000530" "ENSMUSG00000000631" "ENSMUSG00000000817"

Creating enrichGO object

  1. There are ontology options: [“BP”, “MF”, “CC”]. The ont argument can accept either “BP” (Biological Process), “MF” (Molecular Function), and “CC” (Cellular Component) subontologies, or “ALL” for all three.

  2. There are keytypes options: This is the source of the annotation (gene ids). The options vary for each annotation. In the example of org.Mm.eg.db, the options are:

keytypes(org.Mm.eg.db)
##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
##  [6] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"    
## [11] "GENETYPE"     "GO"           "GOALL"        "IPI"          "MGI"         
## [16] "ONTOLOGY"     "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
## [21] "PROSITE"      "REFSEQ"       "SYMBOL"       "UNIPROT"

3.The different organisms with annotation databases available to use with for the OrgDb argument can be found below: Org databases

#go_enrich <- enrichGO(gene = sig genes,
                      #universe = background genes,
                      #OrgDb = organism Db, 
                      #keyType = 'ENSEMBL',
                      #readable = T,
                      #ont = "BP",
                      #pvalueCutoff = 0.05, 
                      #qvalueCutoff = 0.05)

This next part will take a few minutes to run!

#saveRDS(go_enrich, file = "enrichGO_cd8.rds")

go_enrich <- readRDS("enrichGO_cd8.rds")
## Output results from GO analysis to a table
cluster_summary <- data.frame(go_enrich)

write.csv(cluster_summary, "clusterProfiler_TCF7ko.csv")

head(cluster_summary)
##                    ID                                  Description GeneRatio
## GO:0007186 GO:0007186 G protein-coupled receptor signaling pathway    57/661
## GO:0002250 GO:0002250                     adaptive immune response    57/661
## GO:0045785 GO:0045785         positive regulation of cell adhesion    55/661
## GO:0050865 GO:0050865                regulation of cell activation    67/661
## GO:0032943 GO:0032943               mononuclear cell proliferation    43/661
## GO:0032944 GO:0032944 regulation of mononuclear cell proliferation    37/661
##              BgRatio       pvalue     p.adjust       qvalue
## GO:0007186 341/12632 3.388225e-15 1.464052e-11 1.161270e-11
## GO:0002250 363/12632 5.566961e-14 1.202742e-10 9.540014e-11
## GO:0045785 349/12632 1.354755e-13 1.928068e-10 1.529322e-10
## GO:0050865 485/12632 1.784835e-13 1.928068e-10 1.529322e-10
## GO:0032943 248/12632 2.700771e-12 2.262129e-09 1.794295e-09
## GO:0032944 192/12632 3.659525e-12 2.262129e-09 1.794295e-09
##                                                                                                                                                                                                                                                                                                                                                                                                                  geneID
## GO:0007186                                                                   Ccl3/Tulp3/Ccl4/Plek/Adcy1/Pld2/Syk/Iqgap2/App/Rgs11/Atp2b4/Rgs16/Xcl1/Car2/Gnb4/P2ry1/Adgrl2/Adgrb2/Rgs12/Trpm1/Syp/Gpr83/Pde4a/Vipr1/Acpp/Cish/Mgll/Lpar6/Dgkh/Ccl5/Ptger2/Ccr7/Gpr34/Lrp1/Ramp3/Grk3/Gng2/Ccrl2/Gpr68/Baiap3/Cxcr5/Ccr4/Entpd1/F2r/Cxcr6/Mrgpre/Ccr2/Olfr524/Cx3cr1/Rnf157/Pde5a/Adgrg3/Adgrg5/Tcp11/Fzd10/Kctd12/Pde2a
## GO:0002250                                                Tbx21/Cd79a/Fcer2a/Mef2c/Sh2d1a/Cd247/Slamf6/Gzmb/Gata3/Pdcd1lg2/Lyst/Syk/Cd4/Cd74/Cd7/Ctla4/Il18rap/Il18r1/Xcl1/Cr2/Card9/Il6ra/Lef1/P2rx7/Lag3/Cd19/Swap70/Cd40lg/Crtam/Nedd4/Eomes/Cd226/Fut7/H2-Aa/Unc93b1/Prf1/H2-DMa/Ccr7/Prdm1/Slamf7/Cd160/Susd4/Stx11/Foxp3/Fgl2/Cd79b/Serpina3g/Serpinb9/Cd24a/Ccr2/Themis/Cx3cr1/Gzmm/Ifng/H2-Eb1/Tnfrsf13c/H2-Ab1
## GO:0045785                                                                            Cdh1/Pdgfb/Tnfsf14/Fbln1/Slc4a1/Csf1/Itga2/Gata3/Pdcd1lg2/Cyth3/Il12rb2/Myb/Pld2/Gcnt2/Syk/Il6st/Hes1/Cd4/Cd74/Ptprj/Nrp1/Xcl1/Il2ra/Il6ra/Vcam1/Lef1/Nr4a3/Epha1/Ret/Cd40lg/Il15/Pik3r2/Smad3/Chst2/Ccl5/Hpse/Fut7/H2-Aa/H2-DMa/Sirpa/Ccr7/Cd160/Zfhx3/Adk/Foxp3/Prkce/Cd24a/Ndnf/Ccr2/Sox12/Ifng/H2-Eb1/Zbtb16/Tnfrsf13c/H2-Ab1
## GO:0050865 Pdgfb/Myo18a/Tbx21/Hmox1/Mef2c/Tnfsf14/Slc4a1/Blk/Gata3/Pdcd1lg2/Il12rb2/Lyst/Myb/Plek/Id2/Pld2/Syk/Il6st/Dlg5/Hes1/Cd200/Cdkn1a/Cd4/Cd74/Ptprj/Ctla4/Xcl1/Il2ra/Il6ra/Vcam1/Lef1/Nr4a3/Lag3/Cd22/Tyrobp/Cd19/Cd40lg/Fgfr1/Il15/Ubash3b/Crtam/Tirap/Ldlr/Cd226/Ccl5/H2-Aa/H2-DMa/Sirpa/Ccr7/Prdm1/Slamf7/Cd160/Adk/Foxp3/Fgl2/Zc3h12d/Ctla2a/Cd24a/Ccr2/Tlr6/Sox12/Pde5a/Ifng/H2-Eb1/Zbtb16/Tnfrsf13c/H2-Ab1
## GO:0032943                                                                                                                                                  Cd79a/Mef2c/Slc4a1/Blk/Csf1/Slamf6/Pdcd1lg2/Il12rb2/Syk/Il6st/Dlg5/Hes1/Cdkn1a/Cd4/Cd74/Csf1r/Ctla4/Xcl1/Cr2/Il2ra/Vcam1/Lef1/P2rx7/Cd22/Tyrobp/Cd19/Cd40lg/Il15/Crtam/Tirap/Ccl5/H2-Aa/Ccr7/Prdm1/Adk/Foxp3/Zc3h12d/Cd24a/Ccr2/Pde5a/Ifng/Tnfrsf13c/H2-Ab1
## GO:0032944                                                                                                                                                                                   Mef2c/Slc4a1/Blk/Csf1/Pdcd1lg2/Il12rb2/Syk/Il6st/Dlg5/Hes1/Cdkn1a/Cd4/Cd74/Csf1r/Ctla4/Xcl1/Il2ra/Vcam1/Cd22/Tyrobp/Cd40lg/Il15/Crtam/Tirap/Ccl5/H2-Aa/Ccr7/Prdm1/Adk/Foxp3/Zc3h12d/Cd24a/Ccr2/Pde5a/Ifng/Tnfrsf13c/H2-Ab1
##            Count
## GO:0007186    57
## GO:0002250    57
## GO:0045785    55
## GO:0050865    67
## GO:0032943    43
## GO:0032944    37

Visualizing clusterProfiler results

clusterProfiler has a variety of options for viewing the over-represented GO terms.

Bar plot is the most widely used method to visualize enriched terms. It depicts the enrichment scores (e.g. p values) and gene count or ratio as bar height and color. Users can specify the number of terms (most significant) or selected terms to display via the showCategory parameter.

barplot <- barplot(go_enrich, 
                   showCategory = 10,
                   title = "GO Biological Pathways",
                   font.size = 8)

barplot

Dot plot is similar to bar plot with the capability to encode another score as dot size. Go ahead and convert the figure into a dotplot to display the top 10 GO terms by gene ratio (# genes related to GO term / total number of sig genes).

?dotplot
## Help on topic 'dotplot' was found in the following packages:
## 
##   Package               Library
##   enrichplot            /users/m/m/mmg23200/R/x86_64-pc-linux-gnu-library/4.2
##   clusterProfiler       /users/m/m/mmg23200/R/x86_64-pc-linux-gnu-library/4.2
##   graphics              /usr/local/lib/R/library
##   lattice               /usr/local/lib/R/library
## 
## 
## Using the first match ...
#from enrichplot library
dot <- dotplot(go_enrich, 
               showCategory = 10, 
               font.size = 8)
dot

We can also visualize enriched GO terms as a directed acyclic graph:

pathway <- goplot(go_enrich, showCategory = 10)
pathway

Gene Concept Network

Both the barplot() and dotplot() only displayed most significant terms, but some users may want to know which genes are involved in these significant terms.

In order to consider the potentially biological complexities in which a gene may belong to multiple categories, we can use the cnetplot() function. The cnetplot() depicts the linkages of genes and biological concepts (e.g. GO terms or KEGG pathways) as a network.

cnetplot(go_enrich, 
         showCategory = 5, 
         categorySize="pvalue", 
         foldChange=background_gene_list, 
         vertex.label.font=6, 
         cex_label_gene = 0.5, 
         cex_label_category= 0.8)

cnetplot(go_enrich, 
         showCategory = 5, 
         categorySize="pvalue", 
         foldChange=background_gene_list, 
         vertex.label.font=6, 
         cex_label_category = 0.5,
         node_label="category")

Select which labels to be displayed. Options include ‘category’, ‘gene’, ‘all’(the default) and ‘none’.

library(ggplot2)
plot_color <- cnetplot(go_enrich, 
                       showCategory = 5, 
                       categorySize="pvalue", 
                       foldChange=background_gene_list, 
                       vertex.label.font=6, 
                       cex_label_gene = 0.5,
                       cex_label_category= 0.8) + 
  scale_color_gradient2(name='fold change', low='darkgreen', high='firebrick')

plot_color

cnetplot(go_enrich, 
         showCategory = 5, 
         categorySize="pvalue", 
         foldChange= background_gene_list, 
         circular = TRUE, 
         colorEdge = TRUE,
         cex_label_gene = 0.35, 
         vertex.label.font=6)

Circular plots can be quite powerful, if the list was smaller, here things get a bit crowded - overall makes it hard to interpret. Moral of the story, just because you can make it doesn’t mean you should!

tree_cd8 <- pairwise_termsim(go_enrich)
treeplot(tree_cd8)

KEGG Pathway Enrichment

For KEGG pathway enrichment using the gseKEGG() function, we need to convert id types. We can use the bitr function for this (included in clusterProfiler). It is normal for this call to produce some messages / warnings.

# Convert gene IDs for enrichKEGG function
# We will lose some genes here because not all IDs will be converted
ids<-bitr(names(original_gene_list), 
          fromType = "ENSEMBL", 
          toType = "ENTREZID", 
          OrgDb="org.Mm.eg.db") 

# remove duplicate IDS (here I use "ENSEMBL", but it should be whatever was selected as keyType)

dedup_ids = ids[!duplicated(ids[c("ENSEMBL")]),]
# Create a new dataframe df2 which has only the genes which were successfully mapped using the bitr function above
df2 = resdata[resdata$input_id %in% dedup_ids$ENSEMBL,]
# Create a new column in df2 with the corresponding ENTREZ IDs
df2 <- merge(resdata, dedup_ids, by.x = "input_id", by.y = "ENSEMBL")

str(df2)
## 'data.frame':    14850 obs. of  15 variables:
##  $ input_id        : chr  "ENSMUSG00000000001" "ENSMUSG00000000028" "ENSMUSG00000000037" "ENSMUSG00000000056" ...
##  $ baseMean        : num  3904.96 66.3 4.94 833.78 25622.4 ...
##  $ log2FoldChange  : num  0.076 0.302 0.102 -0.226 0.36 ...
##  $ lfcSE           : num  0.0763 0.3278 1.2046 0.1112 0.2128 ...
##  $ stat            : num  0.9957 0.9204 0.0843 -2.0294 1.6914 ...
##  $ pvalue          : num  0.3194 0.3574 0.9328 0.0424 0.0908 ...
##  $ padj            : num  0.593 0.629 NA 0.149 0.26 ...
##  $ WT_CD8_rep1     : num  3.60e+03 5.89e+01 9.99e-01 9.47e+02 1.80e+04 ...
##  $ WT_CD8_rep2     : num  4013.92 53.38 9.53 850.25 25351.27 ...
##  $ WT_CD8_rep3     : num  3795.82 64.92 5.41 899.18 24233.5 ...
##  $ TCF1_KO_CD8_rep1: num  4075.6 78.88 4.53 737.14 32210.41 ...
##  $ TCF1_KO_CD8_rep2: num  3943.78 56.67 5.06 767.1 27000.26 ...
##  $ TCF1_KO_CD8_rep3: num  3998.7 85 4.1 801.8 26937.4 ...
##  $ symbol          : chr  "Gnai3" "Cdc45" "Scml2" "Narf" ...
##  $ ENTREZID        : chr  "14679" "12544" "107815" "67608" ...
# Create a vector of the gene universe
kegg_gene_list <- df2$log2FoldChange
# Name vector with ENTREZ ids
names(kegg_gene_list) <- df2$ENTREZID
head(kegg_gene_list)
##       14679       12544      107815       67608       23849       29871 
##  0.07596003  0.30166718  0.10152331 -0.22561867  0.35991087 -0.14503206
# omit any NA values 
kegg_gene_list<-na.omit(kegg_gene_list)
# sort the list in decreasing order (required for clusterProfiler)
kegg_gene_list = sort(kegg_gene_list, decreasing = TRUE)

We just created the background list, similar to what we did above. Now, we will move on to create the significant list.

# Extract significant results from df2
kegg_sig_genes_df = subset(df2, padj < 0.05)
# From significant results, we want to filter on log2fold change
kegg_sig_genes <- kegg_sig_genes_df$log2FoldChange
# Name the vector with the ENTREZID
names(kegg_sig_genes) <- kegg_sig_genes_df$ENTREZID
head(kegg_sig_genes)
##      74204      14673      18618      12550      67027      54204 
## -0.3414865 -0.3829829 -1.2653917  4.2399463  0.3273819 -0.2018348
# omit NA values
kegg_sig_genes <- na.omit(kegg_sig_genes)
# filter on log2fold change 
kegg_sig_genes <- names(kegg_sig_genes)[abs(kegg_sig_genes) > 1]

BTW, the clusterProfiler package provides search_kegg_organism() function to help searching supported organisms.

sal_ent <- search_kegg_organism('Salmonella enterica', by='scientific_name')
dim(sal_ent)
## [1] 50  3
#sal_ent

Creating an enrichKEGG object

  • organism KEGG Organism Code: The full list is here: https://www.genome.jp/kegg/catalog/org_list.html (need the 3 letter code). I define this as kegg_organism first, because it is used again below when making the pathview plots.

  • keyType one of ‘kegg’, ‘ncbi-geneid’, ‘ncib-proteinid’ or ‘uniprot’.

head(kk)
##                ID
## mmu04060 mmu04060
## mmu04061 mmu04061
## mmu04640 mmu04640
## mmu05321 mmu05321
## mmu05323 mmu05323
## mmu05200 mmu05200
##                                                                                         Description
## mmu04060                        Cytokine-cytokine receptor interaction - Mus musculus (house mouse)
## mmu04061 Viral protein interaction with cytokine and cytokine receptor - Mus musculus (house mouse)
## mmu04640                                    Hematopoietic cell lineage - Mus musculus (house mouse)
## mmu05321                                    Inflammatory bowel disease - Mus musculus (house mouse)
## mmu05323                                          Rheumatoid arthritis - Mus musculus (house mouse)
## mmu05200                                            Pathways in cancer - Mus musculus (house mouse)
##          GeneRatio  BgRatio       pvalue     p.adjust       qvalue
## mmu04060    30/320 129/5552 1.965334e-11 4.451543e-09 3.867418e-09
## mmu04061    18/320  47/5552 3.049002e-11 4.451543e-09 3.867418e-09
## mmu04640    18/320  67/5552 2.161501e-08 2.103861e-06 1.827796e-06
## mmu05321    14/320  46/5552 1.502688e-07 1.096962e-05 9.530204e-06
## mmu05323    12/320  56/5552 6.128943e-05 3.181992e-03 2.764456e-03
## mmu05200    40/320 372/5552 7.111370e-05 3.181992e-03 2.764456e-03
##                                                                                                                                                                                                                                                      geneID
## mmu04060                                                               11482/14103/20302/50930/12977/16162/20303/16195/12504/12978/16174/16182/16963/16184/16194/21949/21947/16168/16154/20304/12775/12145/100532/12773/80901/12772/13051/15978/72049/16188
## mmu04061                                                                                                                                        20302/50930/12977/20303/16195/12978/16174/16182/16963/16184/16194/16154/20304/12775/12145/12773/12772/13051
## mmu04640                                                                                                                                        14128/12977/16398/12504/12978/12482/12516/12902/16184/16194/12483/12478/14960/14998/12484/14969/16188/14961
## mmu05321                                                                                                                                                                57765/14462/16162/16174/16182/24088/17127/14960/14998/20371/15978/17132/14969/14961
## mmu05323                                                                                                                                                                            20302/12977/11975/12477/24088/16168/20304/14960/14998/15978/14969/14961
## mmu05200 12550/18591/14103/15368/16001/17873/16398/16162/432530/18806/69635/16195/15205/12575/18712/58801/12978/21416/226519/16184/13555/14696/16194/16842/19713/14182/16168/18709/13143/17127/67168/19217/18131/14702/14062/15978/21415/235320/16188/93897
##          Count
## mmu04060    30
## mmu04061    18
## mmu04640    18
## mmu05321    14
## mmu05323    12
## mmu05200    40
saveRDS(kk, file = "enrichKEGG_cd8.rds")

kk <- readRDS("enrichKEGG_cd8.rds")
#saveRDS(gseaKEGG_results, file = "gseaKEGG_results.rds")

gseaKEGG_results <- readRDS("gseaKEGG_results.rds")

Visualize enriched KEGG pathways

To view the KEGG pathway, users can use the browseKEGG function, which will open a web browser and highlight enriched genes.

#browseKEGG(kk, 'mmu05202')
library("pathview")
mmu05152 <- pathview(gene.data  = kegg_gene_list,
                     pathway.id = "mmu05152",
                     species    = "mmu",
                     limit      = list(gene=max(abs(kegg_gene_list)), cpd=1))

knitr::include_graphics("mmu05152.pathview.png")

Session info

sessionInfo()
## R version 4.2.2 (2022-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] ggplot2_3.5.0         knitr_1.46            enrichplot_1.18.4    
##  [4] pathview_1.38.0       ggnewscale_0.4.10     DOSE_3.24.2          
##  [7] org.Mm.eg.db_3.16.0   AnnotationDbi_1.60.2  IRanges_2.32.0       
## [10] S4Vectors_0.36.2      Biobase_2.58.0        BiocGenerics_0.44.0  
## [13] clusterProfiler_4.6.2
## 
## loaded via a namespace (and not attached):
##   [1] fgsea_1.24.0           colorspace_2.1-0       ggtree_3.6.2          
##   [4] gson_0.1.0             qvalue_2.30.0          XVector_0.38.0        
##   [7] fs_1.6.3               aplot_0.2.2            rstudioapi_0.16.0     
##  [10] farver_2.1.1           graphlayouts_1.1.1     ggrepel_0.9.5         
##  [13] bit64_4.0.5            fansi_1.0.6            scatterpie_0.2.2      
##  [16] codetools_0.2-18       splines_4.2.2          cachem_1.0.8          
##  [19] GOSemSim_2.24.0        polyclip_1.10-6        jsonlite_1.8.8        
##  [22] GO.db_3.16.0           png_0.1-8              graph_1.76.0          
##  [25] ggforce_0.4.2          compiler_4.2.2         httr_1.4.7            
##  [28] Matrix_1.5-1           fastmap_1.1.1          lazyeval_0.2.2        
##  [31] cli_3.6.2              tweenr_2.0.3           htmltools_0.5.8.1     
##  [34] tools_4.2.2            igraph_2.0.3           gtable_0.3.4          
##  [37] glue_1.7.0             GenomeInfoDbData_1.2.9 reshape2_1.4.4        
##  [40] dplyr_1.1.4            fastmatch_1.1-4        Rcpp_1.0.12           
##  [43] jquerylib_0.1.4        vctrs_0.6.5            Biostrings_2.66.0     
##  [46] ape_5.7-1              nlme_3.1-160           ggraph_2.2.1          
##  [49] xfun_0.43              stringr_1.5.1          lifecycle_1.0.4       
##  [52] XML_3.99-0.16.1        org.Hs.eg.db_3.16.0    zlibbioc_1.44.0       
##  [55] MASS_7.3-58.1          scales_1.3.0           tidygraph_1.3.1       
##  [58] parallel_4.2.2         KEGGgraph_1.58.3       RColorBrewer_1.1-3    
##  [61] yaml_2.3.8             memoise_2.0.1          gridExtra_2.3         
##  [64] downloader_0.4         ggfun_0.1.4            HDO.db_0.99.1         
##  [67] yulab.utils_0.1.4      sass_0.4.9             stringi_1.8.3         
##  [70] RSQLite_2.3.6          highr_0.10             tidytree_0.4.6        
##  [73] BiocParallel_1.32.6    GenomeInfoDb_1.34.9    rlang_1.1.3           
##  [76] pkgconfig_2.0.3        bitops_1.0-7           evaluate_0.23         
##  [79] lattice_0.20-45        purrr_1.0.2            labeling_0.4.3        
##  [82] treeio_1.22.0          patchwork_1.2.0        cowplot_1.1.3         
##  [85] shadowtext_0.1.3       bit_4.0.5              tidyselect_1.2.1      
##  [88] plyr_1.8.9             magrittr_2.0.3         R6_2.5.1              
##  [91] generics_0.1.3         DBI_1.2.2              pillar_1.9.0          
##  [94] withr_3.0.0            KEGGREST_1.38.0        RCurl_1.98-1.14       
##  [97] tibble_3.2.1           crayon_1.5.2           utf8_1.2.4            
## [100] rmarkdown_2.26         viridis_0.6.5          grid_4.2.2            
## [103] data.table_1.15.4      Rgraphviz_2.42.0       blob_1.2.4            
## [106] digest_0.6.35          tidyr_1.3.1            gridGraphics_0.5-1    
## [109] munsell_0.5.1          viridisLite_0.4.2      ggplotify_0.1.2       
## [112] bslib_0.7.0